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Abstract 

Remodelling is denned as an evolution of microstructure or variations 
in the configuration of the underlying manifold. The manner in which a 
biological tissue and its subsystems remodel their structure is treated in 
a continuum mechanical setting. While some examples of remodelling are 
conveniently modelled as evolution of the reference configuration (Case 
I), others are more suited to an internal variable description (Case II). 
In this paper we explore the applicability of stationary energy states to 
remodelled systems. A variational treatment is introduced by assuming 
that stationary energy states are attained by changes in microstructure 
via one of the two mechanisms — Cases I and II. An example is presented 
to illustrate each case. The example illustrating Case II is further studied 
in the context of the thermodynamic dissipation inequality. 

1 Introduction and background 

The development of a biological tissue and its subsystems consists of the distinct 
processes of morphogenesis, growth and remodelling, a classification suggested by 
Taber (1995). For preciseness of mathematical formulation we have previously 
defined and treated growth as consisting of only addition or depletion of mass 
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through processes of transport and reaction, possibly coupled with mechanics 
(Garikipati et al., 2004). We define remodelling as micro structural changes 
within the biological structure at constant mass. While remodelling and growth 
occur simultaneously and in a coupled fashion in biological tissues, they can be 
treated as separate processes for modelling purposes. Furthermore, in certain 
situations, addition and depiction remain in balance, maintaining constant mass. 
This is referred to as homeostasis, during which remodelling can occur. 

There is, in fact, experimental evidence for a strict separation between re- 
modelling and growth in soft tissue: (i) Stopak and Harris (1982) described the 
orientation of collagen fibrils due to the forces exerted on them by the fibroblasts 
(tendon cells) in a collagen gel. Growth processes of resorption of existing fibrils 
and production of new ones with the preferred orientation were not reported 
in their paper, suggesting that it was the mechanical action of fibroblasts alone 
that resulted in fibril orientation. In this example the fibril orientation can be 
viewed as the microstructural quantity undergoing an evolution in the absence 
of growth. Fibril reorientation driven by stress and occurring independently of 
growth was attained in our laboratory: A collagen gel was formed in a dish with 
polyethylene supports fixed to the base. When fibroblasts were added to the gel 
they exerted traction, thereby aligning the collagen fibrils. The orientation that 
was obtained corresponded to the axis defined by the supports, (ii) Another 
instance of microstructural change in collagen fibrils is their longitudinal and 
lateral fusion. Birk et al. (1995, 1997) reported an abrupt change in the length 
of embryonic chicken tendons from ps 40 /im to « 120 /im, on the 17 th day 
after fertilization. In this case, while growth obviously was taking place in the 
embryonic tendon, micrographs verified the abrupt change to be due to longi- 
tudinal (end-to-end) fusion of smaller fibrils. This time scale, over which fusion 
took place, is much smaller than the time scale of growth. We therefore propose 
that, in a mathematical treatment, it is appropriate to view this remodelling 
process as taking place in the absence of growth. 

In contrast to these examples in soft tissue stands the case of bone, wherein 
the macroscopic process that is usually labelled "remodelling" takes place by 
resorption of existing collagen fibrils and production of new fibrils with a pre- 
ferred orientation. In this case, therefore, the finer scale processes do fit our 
definition of growth. For the sake of conceptual and mathematical clarity we 
will ignore all processes that require growth at any scales in this paper, and 
focus upon a continuum mechanical treatment of remodelling. 

The literature in biomechanics and mathematical biology has a number of 
papers concerning "remodelling" (Cowin and Hegedus, 1976; Taber, 1995; Har- 
rigan and Hamilton, 1993; Scliktar et al., 2000; Taber and Humphrey, 2001; 
Ambrosi and Mollica, 2002). Almost universally, however, these papers treat 
mass and density changes and the mechanics — characterized by internal stress — 
that is associated with them. Therefore, by our definition, they describe growth. 
An exception is Humphrey and Rajagopal (2002), in which the notion of a nat- 
ural configuration is introduced. It bears some similarities with the treatment 
in Section 2 of this paper. 

Another treatment, by Driessen et al. (2003), is specific to the evolution of 
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fiber orientation within a tissue. However, there are at least two fundamental 
shortcomings in their model: (i) It is unable to distinguish between cubic or- 
thotropy and isotropy, and (ii) it predicts that in a tissue that is undcformcd 
from its reference state, the fiber orientations evolve until a tensorial variable 
representing their distribution reaches a particular value. We discuss their model 
and present a critique of both these features in Section 4.4. 

The following is the organization of this paper: The mathematical treat- 
ment for cases that are best described by smooth configurational changes on 
the underlying manifold (Case I remodelling) is presented in Section 2. A one- 
dimensional example is presented in Section 3 to illustrate this formulation. 
The local reorientation of collagen fibrils, including our experiments and the 
treatment via internal variables is discussed in Section 4. Thermodynamic dis- 
sipation is discussed in Section 5. The paper concludes with a brief discussion 
in Section 6. 

2 Case I remodelling: Microstructural changes 
that alter the reference configuration 

Figure 1 depicts the kinematics associated with Case I remodelling. Material 
particles arc labelled by X £ I 3 in the reference configuration, which is de- 
noted by fio C M. 3 . The material microstructure undergoes changes that can 
be described by a point-to-point vector map, x- x [0,T] i— > R 3 , defined as 
x{X,t) = X + n(X,t). It carries the microstructure from the reference con- 
figuration to a remodelled configuration, 12*, in which material points will also 
be labelled as X* = X- Assuming t) to be smooth in X, its tangent map 

is K(X,t) = dx/dX, leading to K = 1 + dn/dX. In Case I remodelling, 
therefore, K denotes a compatible change in configuration. 1 

Distinct from x is the point-to-point vector map (p* : O* x [Q,T] i— > K 3 . 
It carries material points from 12* to the spatial configuration 12, and is the 
deformation relative to 12*. The placement of material points in 12 is therefore 
x = (p*(X*,t). The displacement, u* , satisfies cp*(X*,t) = X* + u*(X* ,t). 
The classical deformation gradient is F* := dip* /dX* = 1 + du* /dX* . In 
this initial treatment we do not consider any further decompositions of F . 
The overall motion of a point is tp(X,t) = n(X,t) + u*(X*,t) o x{X,t), and 
the corresponding tangent map is F = dtp/dX. It admits the multiplicative 
decomposition F = F*K. 

Remark 1: In applications, the microstructural motion, Xi will be distin- 
guished from the deformation, y>*, by either physical or mathematical consid- 
erations. For instance, in the example of Section 3, x denotes the uncoiling 
of a long-chain molecule, while f* is the stretching of bonds. A mathematical 
decomposition can be motivated in other applications such as x denoting fine 
scale motion of material particles, while ip* denotes motion on coarser scales. 

lr This configurational change can be further decomposed, multiplicatively, into incompati- 
ble components: K = K X K 2 as in Garikipati et al. (2005). 
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Figure 1 : The kinematics of Case I remodelling and deformation. 
2.1 A variational formulation 

We wish to explore the relevance of stationary energy states to remodelled 
biological systems. For this purpose the microstructural changes of the type 
outlined with examples from biology in Section 1 and discussed more math- 
ematically in this section, are assumed to occur when the Gibbs free energy 
of the system attains a stationary state. Importantly, the stationary state of 
the biological system is not one of equilibrium. The latter class of states is 
dictated by thermodynamic dissipation (see de Groot and Mazur, 1984). The 
stationary energy states, on the other hand, can be identified via standard vari- 
ational arguments. For this purpose we consider the following Gibbs free energy 
functional: 

Ui[u*,k] = j ii*{F*,K,X*)dV* 

- J f*-(u* + K)dV* - J i*-(u* + #c)dA*, (1) 

where ip* = ip* (F* , K, X*) is the Helmholtz free energy density function, de- 
fined per unit volume in the remodelled configuration. Observe that ip* is 
assumed to depend upon K in addition to the usual dependence on F* . Ma- 
terial heterogeneity is allowed, and is represented by the dependence of ip* on 
X* . The body force per unit volume in fi* is /*, and the applied traction per 
unit area of the surface subset, dVl^ C dCl*, is i . We allow for displacement 
boundary conditions, u* = g* . and vanishing microstructural change, k = 
on 9fi* = <9£T\9ilJ\ Since the total motion of a material point is k + u* , the 
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potential energy of the external loads is as expressed by the second and third 
terms in (1). 



2.1.1 Euler-Lagrange equation: Quasistatic balance of linear mo- 
mentum 

The functional (1) yields one set of Euler-Lagrange equations when stationarity 
of 11/ is imposed with respect to variations in u* . We first define u* := u*+e6u* , 
at fixed k, with 5u* £ M 3 and Su* = on cftl* . Then, 

^HjK.k] = 1 Ur(F* £ ,K,X*)dV*- J f*-(u* s + n)dV* 
\i* n* / 

~ J ?.(u* e + K)dA*. 

an* 

Differentiating under the integrals, applying the chain rule, integrating by parts, 
and imposing stationarity via (dII//d£) e= o = gives 

- / Div*[^]-5u*dV* - J f*-5u*dV* 

Q* n* 

+ J (S^mA -5u*dA* - J ?-5u*dA* = 0. (2) 

an* an* 

Introducing P* = difj*/dF*, the arbitrariness of Su* £ M 3 and the standard 
localization argument yield the Euler-Lagrange equation 

Div*P* + /* = Oinfr, (3) 

the following boundary condition and constitutive relation 

P*N*=i* on dn* t> P* = ^- (4) 

Of these ccmations, (3) is recognized as the quasistatic balance of linear momen- 
tum in Q* , and (4)i as the corresponding traction boundary condition. Observe 
that, as defined in (4)2, P* is the first Piola-Kirchhoff stress, which is conjugate 
to F* with tp* as the relevant strain energy density function. 

2.1.2 Euler-Lagrange equation: Stationarity with respect to microstruc- 
tural change 

The class of variations now considered is given by 

K e := k + sSk, at fixed it*, with^K = on9f2*. (5) 
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Proceeding as in Section 2.1.1: 



U t [u*,k £ ] = A ^J$*(F* s ,K e ,X* E )dV*- J f*-(u* + K £ )dV* 

J i*-(u* + K e )dA*, 



a* 

d 
dl 



where the following relations hold: 



dS k 



f; = f £ k; 1 , x; = x + K e , k £ = i + ok/ox + s— . (6) 

OX 

Standard, if lengthy, manipulations (see Appendix A) then lead to the fol- 
lowing set of equations governing the stationary energy state in which the con- 
figurational variables take on the values k — k s and K = K s . 



-Div* (^*1-F* T P* + S*) + ^ =0 in Q* (7) 

(ip*l — F* T P* + S*^ N* = on dCl*, (8) 

where S* = -^K T in fi* (9) 

T 

Observe that the Eshclby stress ip*l ~ F* P* makes its appearance. Hereafter, 
it will be denoted by £. The quantity X* is a thermodynamic driving force 
defined in (9) as the change in Hclmholtz free energy density corresponding to a 
change in the tangent map of the microstructural configuration. It is stress-like 
in its physical dimensions and is a second-order tensor. For this reason we refer 
to it as a non-Eshelbian configurational stress. The boundary condition in (8) 
is simply a restriction on the normal component of this thermodynamic driving 
term. 



2.1.3 Stationarity with simultaneous variation of u* and k 

The procedure followed above is formal in the sense that the independent im- 
position of variations has been assumed: u* = u* + eSu* for 5k = 0, and 
K £ = K + e5k for Su* = 0. Physically, this may not be possible due to the in- 
teraction of deformation and configurational changes. Stationarity of the Gibbs 
free energy under simultaneous variation of its arguments is obtained by requir- 
ing 

— IL I [vl,K e ]\ e=o =0. 

On carrying out this calculation as in Sections 2.1.1, 2.1.2 and Appendix A 
it can be shown that stationarity requires the following condition: 
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(Div*P* + /*) -Su* + (^Div* (£ + E*) + -Snj dV* 

n* 

- J ((P*N* - T) -Su* + ((£ + E) AT*) -6k) dA* = 0, (10) 

an; 

where the notation introduced above for the Eshelby stress and the non-Eshclbian 
configurational stress has been used. Clearly (3-4) and (7-9) ensure satisfaction 
of (10), and are therefore sufficient conditions for stationarity. 

3 An example of Case I remodelling: Configu- 
rational change of a long chain molecule 

In this section we employ an established statistical mechanical example of con- 
figurational changes of long chain molecules in order to illuminate the math- 
ematical formulation of Section 2. In addition to the interest in this example 
from the standpoint of this paper we point out that configurational changes are 
of central importance to the chemical activity of long chain molecules. 

Consider a long-chain molecule, say a protein, that can exist in a highly 
coiled state. Let its contour length be L. In the reference state, Qq, the end- 
to-end lengths of the coiled and the relatively straight domains are in the ratio 
£ : 1 — £, and the end-to-end length of the molecule is ro (Figure 2). An 
entropic elasticity is associated with uncoiling of the molecule. Let k be the 
change in end-to-end length from the reference state due to uncoiling; it is a 
measure of changes in the number of configurations available (entropy). The 
corresponding contribution to the Gibbs free energy function is specified by the 
worm-like chain (WLC) model (Kratky and Porod, 1949). A stiffness, /i, with 
physical dimensions of energy is associated with bond stretching. The increase 
in length due to bond stretching is u* . The square of the ratio of this quantity 
and the length of the molecule, k + ro, determines the energy stored in bonds. 
The enthalpic elasticity is due to this mechanism. The molecule is subjected to 
an externally-applied axial tensile force, T. 

The one-dimensional tangent map of the local configurational change, av- 
eraged over the entire molecule, is K = 1 + k/tq. It carries the molecule to 
its remodelled configuration, f2*. The one-dimensional deformation gradient 
relative to f2* is F* = 1 + u*/(k + ro), also averaged over the molecule. It 
carries the molecule to its current configuration, f2, in which it has undergone 
a configurational change and bond stretching. 

We will examine the response of the molecule to the axial force. In this case, 
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Figure 2: A long chain molecule that can undergo configurational changes and 
bond stretching due to an external axial load. For ease of visualization the 
figure depicts bond stretching only in the uncoiled, straight regions. 



the Gibbs free energy is 



ILi[u*,k] = ^/i 



2 \ k + r 

ksO ( (k + r ) 2 L K + r 



A \ 2L 4(1- (« + r )/£) 
-T(u* + k), (11) 

where the first term on the right hand-side is the Helmholtz free energy from 
elastic stretching. The second term is the entropic contribution from the WLC 
model, after Marko and Siggia (1995), with ks being the Boltzmann constant 
and 9 being the temperature. The persistence length is A and is defined as the 
ratio of bending stiffness to the thermal energy. It is also the distance along the 
molecule's contour over which the correlation between tangent vectors falls to 
e _1 . See Landau and Lifshitz (1951) for the statistical mechanics behind these 
aspects of the model. The third term in (11) is the potential of the external 
force. 

Proceeding as in Section 2 wc first seek stationarity with respect to variations 
in u* 

u" 

For this example, in the absence of a body force (3) reduces to the trivial 
requirement that the axial force be constant along the molecule. Equation (12) 
is the traction boundary condition (4)i for the present case. Solving, we get 

u* s = -( K + r Q ) 2 , (13) 

where superscript (•)" denotes a quantity with respect to which the system is 
in a stationary state. 



de 



—TLi[u*,k] „ =m ,„ , TT2 ~T = 0. (12) 
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Next, considering variations with respect to k and imposing stationarity 
gives 



-^H/ [u*,K s ] 



e=0 




A(l-( K + r )/Lj 2 



(14) 

In this case the differential equation (7) is trivially satisfied since the forces 
corresponding to the stresses £ and X* are constant along the molecule. 

The first term in the second member of (14) arises due to the variation of F* 
with k. It represents the effect of variation in the underlying manifold, due 
to configurational changes. On comparison with (38) and (39), and following 
the derivation in Appendix A it is clear that this term is the reduced version 
of the Eshelby stress for the present one-dimensional setting. The third term of 
the second member of (14) arises due to the variation of the WLC term with k. 
It represents the non-Eshclbian stress reduced to this setting as is confirmed by 
comparison with the development in Appendix A. 

The reduced version of the boundary condition (8) is obtained on substitut- 
ing the stationary solution (13) for u* in Equation (14). (In effect, this is the 
substitution referred to at the end of Appendix A.): 



{K + r ) + —\ — h — 77772-7 )-T = 0. (15) 

M A y L 4 (1 - (k + r )/Ly 4 J 

Solutions to (15), denoted by k s , can be obtained in closed-form since it is a 
cubic equation. In order to illustrate the nature of the solution we have plotted 
the left hand-side of (15) against k in Figure 3 with the numerical values of 
parameters in Table 1. The plot is restricted to the interval < k < L — ro, 
since uncoiling is of interest, and (11) is non-physical for k > L — tq. 



Tabic 1: Parameters used in the Gibbs free energy of the long chain molecule 



P arameter 


Value 


Units 


Notes 




1 









300 


K 




A 


14.5 


nm 


Collagen monomer molecule (Sun et al., 2002) 


ro 


120 


nm 




L 


309 


nm 


Collagen monomer molecule (Sun et al., 2002) 


/'• 


2.951 x 10~ 7 


J 


Stiff bonds 


T 


10 


pN 


Force applied in Sun et al. (2002) 



The variation in free energy DHj = (dn//de) e= o is negative over most of 
the K-interval and passes through zero at k s = 162.642 nm. This is the length 
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Dili (N) 



1x10' 



,-10 




5x10" 



r" 



IxlU 7 UsxIO 7 lTxTcr 7 1.75x10" 



I" 7 



k (nm) 



/.5x 10 



-5xl0~ 



-1x10" 



,-10 



Figure 3: Variation of DUj = (dri//de) e= o with respect to k. 

increase due to uncoiling when the molecule is in the stationary energy state. 
Corresponding to this value is the bond stretch, u* s = 3 x 10~ 9 nm. The 
molecule is very stiff — almost rigid — to bond stretching. Also observe that the 
stress-stretch response determined by the WLC model locks as k s — > L — ro 
from below (k s — > 189 nm); i.e., as the uncoiled length approaches the contour 
length. 

Remark 2: The uncoiling, represented by re, takes place under an axial force 
and its effect on the free energy is treated via a model of entropic elasticity. 
It therefore follows that the entropy of the molecule is reversibly changed by 
application and removal of a force. 

4 Case II remodelling: Microstructural changes 
represented by internal variables 

While microstructural changes always imply an evolution of the reference con- 
figuration, a simpler description may be possible than the general one developed 
in Section 2. Such an instance is illustrated by the example of collagen fibril re- 
orientation under mechanical load: The reorientation of a fibril does modify the 
reference configuration since the underlying microstructure changes, and hence 
the general formulation of Section 2 remains applicable. However, the physics 
of the process admits a simpler model that requires only the introduction of a 
rotation tensor as an internal variable and is discussed below. We begin with 
the biological basis and include a description of our experiments. 

4.1 Collagen fibril alignment in gels by cell traction 

Collagen is the most widely-present protein in the human body. It exists in many 
forms of which, to fix ideas, we will consider type I collagen. It has a fibrillar 
structure and is the main form of collagen in tendons. The fundamental unit 
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is the triple helical collagen molecule that is assembled from three single chain 
molecular strands. The triple helix has been reported to be between 300 to 360 
nm in length and 1.5 nm in diameter. The triple helices further assemble into 
collagen fibrils of varying lengths (the range of 20-140 /xm has been reported) 
and 10-300 nm in diameter. The fibrils assemble into fibers that can range up 
to millimeters in length. The collagen fibrils arc surrounded by a dense network 
of proteoglycan (PG) molecules, which, at one end, associate with the collagen 
fibrils. They form a hydrated gel that contains most of the fluid phase of the 
extracellular matrix. Other extracellular matrix proteins are also found in the 
gel. The interested reader is directed to Alberts et al. (2002, chap. 19) and 
references therein for details. 

Stopak and Harris (1982) have described the reorientation of collagen fibrils 
due to cell traction in gels. In their experiments, cell-bearing tissue explants 
were embedded in collagen gels. Cells were found to migrate outward from 
the explant, and, by applying traction, to induce a preferred alignment of the 
collagen network. Specifically, collagen fibrils were found to align either between 
a pair of explants, or an explant and a fixed, effectively rigid, support. In a three- 
dimensional environment, such as a collagen gel, the fibroblasts themselves align 
along directions of maximum principal tensile stress in the matrix (Balaban 
et al., 2001). In vitro (and probably in vivo also), the stress is often imposed by 
fibroblast traction on the matrix. 2 In other in vitro studies an externally-applied 
traction has resulted in fibroblast alignment with the maximum principal tensile 
stress direction. 

Fibroblasts attain alignment by attaching themselves to the matrix (usually 
the collagen fibrils) by "three-dimensional adhesion points" . Complex signalling 
pathways and chemical cascades are involved in the formation of these adhesion 
points (Cukicrman et al., 2001). The actual attachment to the matrix is me- 
diated by receptors called integrins that pass through the cell membrane and 
are attached to the intracellular actin network at the opposite end from the 
adhesion point. See Geiger et al. (2001); Mitra et al. (2005). A tensile stress 
develops in the actin network, the cell's interinsic contractile apparatus, asso- 
ciated with stretching of the fibroblast between multiple adhesion points in the 
extra-cellular matrix. The stress in the actin network thus enables the fibroblast 
to apply traction to the matrix in a direction now determined by the alignment 
of adhesion points (Balaban et al., 2001). This traction induces a marked re- 
orientation among the fibrils in a network with an initially random orientation 
distribution. The fibrils align with the maximum principal tensile stress di- 
rection in the matrix. There now exists general agreement that these mutual 
interactions between stress in the matrix, fibroblast alignment and stress in the 
actin network are responsible for the collagen fibril reorientation described in 
papers such as Stopak and Harris (1982). 

2 The aligned fibrils stiffen the matrix in their direction, thereby providing greater resistance 
to fibroblast traction. A "positive-feedback loop" is thus created between fibroblast and 
collagen fibril alignment, a phenomenon that is sometimes called "contact guidance" (Barocas 
and Tranquillo, 1997). 
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4.2 Experimental study 

In order to assess the role of fibroblasts in the alignment of collagen fibrils we 
carried out a set of simple collagen gel contraction experiments (that we have 
already alluded to in Section 1). 

Three experimental configurations were considered. In the first, a cell- 
collagen solution was delivered over porous polyethylene posts spaced 6 mm 
apart (Figure 4). Over 12 hours, the cell-seeded gel contracted to form a lin- 
ear structure stretched between the supporting posts. The contraction of the 
cell-collagen gel was observed (at ambient temperature 22°C) between crossed 
polarizing lenses mounted on an inverted microscope for 12-15 hours under 40x 
magnification using an attached digital camera. An observed increase in trans- 
mitted light through the polarizing filters over the recording period indicated 
a significant microstructural reorganization in addition to macroscopic changes 
in gel shape. The gel became birefringent due to alignment of collagen fibrils 
along the axis of the posts. In a second experiment, a nearly identical proce- 
dure was followed but without the addition of the constraining posts to the dish. 
While this gel contracted to half its original diameter overnight, the lack of con- 
straints resulted in an isotropic contraction, no observable alignment of fibrils, 
and therefore no birefringence (data not shown). Finally, a third gel plated 
without cells was found to neither contract nor yield a birefringent species. 
Taking these results together, one can (reasonably) conclude that the alignment 
requires both traction-providing cells (i.e., fibroblasts) as well as fixed displace- 
ment constraints so that the stress field applied by the cells to the matrix is 
anisotropic and induces the reorientation of fibrils. 




100 mm 



Figure 4: Schematic figure of the experimental setup for collagen gel formation 
and fibril alignment by cell traction. 
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4.3 Modelling assumptions leading to Case II 

We make the following modelling assumptions on the alignment of fibrils by cell 
traction in a collagen gel: 

1. We restrict ourselves to situations in which there is a uniform distribution 
of fibroblasts and fibrils in the gel. The typical linear dimension of a 
fibroblast is 25 /im. Therefore, at a macroscopic scale, each point in the 
gel has available the microstructural mechanisms required for collagen 
fibril reorientation. 

2. Since the fibroblasts and fibrils arc uniformly distributed, cell traction 
does not result in translation of the center of mass of the collagen fibril 
network within an infinitesimal neighborhood of a point. Therefore, the 
only effect of fibroblast traction on the fibrils is to change their alignment 
at a continuum point. 

3. A consequence of Assumption 2 is that the fluid surrounding the collagen 
fibrils in the gel does not undergo a translation of its center of mass due 
to fibroblast traction. 

4. There is no loss of contact between the fibrils and the surrounding fluid 
in the gel due to fibril realignment; i.e., local compatibility is maintained. 

Reverting to Equation (1) we let the motion of material points induced by 
fibril rotation be denoted by n. It follows from Assumptions 2 and 3 that the 
work done by forces /* and i* on k vanishes in (1) . 

At the macroscopic scale, we develop an internal variable description of the 
reorientation of the microscopic fibrils. In contrast to Case I remodelling (Sec- 
tion 2) we do not undertake a detailed consideration of the changes in reference 
configuration brought about by variation of this internal variable. The Gibbs 
free energy will depend on the internal variable for fibril reorientation since it 
determines pointwise anisotropy of the extracellular matrix. If the tangent to 
a fibril has an initial orientation given by the unit vector field Mq(X) 6 Oo, 
then M*(X, t) = Q(X, t)Mo(X) is its remodelled orientation at time t, where 
Q(X,t) € SO(3), a rotation tensor, is the internal variable. We write the 
Hclmholtz free energy density with respect to f2o as ip(F,Q,X). Note that ip 
is defined with respect to the reference configuration, changes in which we will 
ignore in Case II remodelling. 

The Gibbs free energy functional is parametrized by the displacement u, 
now defined on ^o, and Q. Reflecting the assumptions made above and the 
arguments following from them, Equation (1) reduces to 



Note that the integrals are defined on Qq. The body force and traction 
vectors are / and t, respectively, defined per unit volume in f2o and per unit 




(16) 
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surface area on dflo. We now consider variations on Q, recalling that 50(3) is 
a Lie group whose Lie algebra consists of real, skew-symmetric matrices. This 
Lie algebra is denoted by Skw(3) for the case wherein Q acts on vectors in M 3 
(see any standard textbook on manifold analysis, for instance, Choquct-Bruhat 
et al. (1982)). The admissible variations on Q are therefore of the form 



Q e = Q + SQ, where 5Q = WQVW G Skw(3). 
This leads to the following variation on II//: 



(17) 



-U H [u,Q £ ]\ e=Q 



— | / ^(F,Q E ,X)dV 



f-udV 



t-udA 



:=0 



(18) 



dVtt 



Imposing stationarity, invoking the localization argument and applying (17) this 
reduces in a straightforward manner to 



dtp 
dQ 



WQ = 0, V W G Skw(3), 



(19) 



a relation also obtained by Vianello (1996), albeit without beginning from the 
integral form. As a final step this result for a stationary state can be rewritten 

as 

dtp T 



dQ 



Q 1 



Q=Q S 



W = 0, V^eSkw(3). 



(20) 



4.3.1 Example: Collagen fibril reorientation in a gel using a contin- 
uum strain energy function 

The WLC model (11) is extended to a continuum strain energy function for the 
system consisting of the gel, collagen fibrils and fibroblasts. This is achieved by 
introduction of the fibril number density, N, and addition of a repulsive term 
to enforce a vanishing stress at unit stretch and a bulk compressibility term for 
the gel: 



*KF,Q,X) 



Nk B (2r z L 



AA \ L 1 - r/L 
AA \L 4r (l - r f) 2 4r c ' &{ ' 



repulsive term 

-?( J ~ 2/ '- 1 ) + 2 7l! E. (21) 

P 



bulk compressibility 
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The elastic stretch in the direction of the remodelled fibril is A, and E = 
i(_F T F — 1) is the Lagrange strain. The factors 7 and (3 control bulk compress- 
ibility. The end-to-end length is given by 



Jrp?, A = VM*-F T FM*. 



(22) 



Note that the end-to-end length is modified by the macroscopic deformation 
through (22). This is in contrast with the explicit inclusion of the change in 
molecular configuration, k, in (11). The change in configuration considered in 
the present section is at the more macroscopic scale of fibril reorientation. 

Consider a cylindrically-shaped collagen gel in which the collagen fibrils are 
all radially-oriented at time t = 0. Adopting a cylindrical coordinate system, 
{e^, e ai ez}, we have Mq = e R . The gel is loaded by imposing a deformation 
gradient: F = A^e^, (g) e R ■ + X a e a ®e a + \z&z ® e^, such that \z > 1, A_r = X a 
and Xz > X R . We will assume that Case II remodelling occurs according to 
the cell traction mechanisms discussed above, and that the fibrils undergo local 
reorientation with M* — > ez as t — > 00. The rotation tensor corresponding to 
this stationary state is Q s = — en ® ez + e a <E> e a + ez (E> e^. We aim to verify 
that Q s satisfies (20). 

From (21) and (22) we have, after some manipulations, 



Q=Q S 



W 



Nk B 9r 
4A 



4r 

T 



{l-r/Lf 



4r| 
rL 



■F T F(Q S M Q 



'•() 



r(l-r /L) 2 
g)Q s M ): W 



r 



(23) 



Using the tensor product expansions established above for F and Q s , and Mq 
en, this reduces to 



Q=Q a 



W 



Nk B 9r 
4A 



4r 

T 



{l-r/L)' 



-1 - 



'0 



rL r(l - r /L) 2 
X 2 z e z ®e z : W. 



r 



(24) 



Since (24) involves the scalar product of a symmetric tensor and W € 
Skw(3), we have 



Q=Q 3 



W = 0, VWeSkw(3), 

Q s = -e R g) e z + e a ® e a + e z ® e R . (25) 



This implies that a stationary energy state is achieved when the fibrils with 
initial radial orientation undergo reorientation along the direction of maximum 
principal tensile stretch. 
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Remark 3: Even though the discussion in Sections 4.1 and 4.2 refers to fibril 
alignment with the maximum principal tensile stress direction, it is entirely 
equivalent to consider alignment with the maximum principal stretch direction 
in the stationary energy state. This is on account of the work of Vianello (1996) 
who showed that the strain energy of an anisotropic solid is at a minimum when 
the principal stress and stretch directions coincide. This result is used to justify 
an evolution law for fibril orientation (32) below. 

4.4 Notes on a recent model of fiber orientation 

To conclude this section, we briefly discuss a recent theory of remodelling, de- 
veloped by Driessen et al. (2003), for the evolution of fiber orientation within a 
tissue. Their theory is quite distinct from the development in Sections 4.1-4.3, 
and has at least two fundamental shortcomings: 

(i) In order to represent the distribution of fiber orientation, these authors 
work with a fiber orientation tensor, defined as So = (eo, eo) := X)i=i V^o® 
e . There are N distinct orientations, each specified by a unit vector, e l , 
and with a probability distribution ipi, such that y^--, ipi = 1. The ten- 
sor, So, however fails to distinguish between cubic orthotropy (N = 3, 
ipi = 1/3, e -e = the Kronecker-delta) and isotropy (ipi — l/N, the 
directions e being distributed uniformly over the unit sphere, in the limit 
N — > oo). A straightforward calculation shows that So = 5I hi both 
cases. This fundamental failing suggests that higher-order statistics must 
be included, for instance, moments of the fiber distribution. 

(ii) The second limitation is related to the evolution law prescribed for the 
fiber orientation tensor in the current placement, S = (1/ A 2 )FSoF T , 
where A 2 is the mean square of the stretches of the fibers: 

N 

i=i 

v v 
The evolution law prescribed is S +2(D : S)S = -^(A — S), where S is the 

Lie derivative of S, D is the rate of deformation tensor, t is the relaxation 

time, and A is written as a function of the finger tensor: A = B v /tr(B u ) 1 

with v being a real exponent. The physical implication of such a law is that 

S continues to evolve, driven by the strain, until S = B v /tv(B y ): For a 

point with B = 1, i.e., unchanged from its reference placement, this means 

that any initial fiber orientation tensor, ■S'o(O) = 5(0) (e.g., So(0) — 

S(0) = e ® e ) must evolve until So{t) = S(t) = 1. We know of no 

physiological instances in which this happens. In fact, this conclusion does 

not correspond with common experience. In our experiments described 

in Section 4.2 we have also confirmed that the fiber orientations in a 

uniaxially-aligned collagen gel remain in this state if the gel does not 
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deform, and do not evolve to S(t) = 1 as Driessen and co-workers' rule 
suggests. 

5 Thermodynamic dissipation associated with 
reorienting fibrils 

We return to our formulation of collagen fibril rc-oricntation via internal vari- 
ables for the discussion on dissipation. Denoting the internal energy of the 
system (gell, collagen fibrils and fibroblasts) by e, and the energy lost against 
viscous resistance by S> v , the First Law of thermodynamics gives 

e = P: F-% - V-q, (26) 

where P is the first Piola-Kirchhoff stress defined on f^o, and q is the heat flux. 
The energy lost against viscous resistance as the fibrils rotate relative to the 
viscous gel is in the form of heat, therefore — & v is heat lost from the system. 
To be more precise we assume a linear viscosity and write &! v = iju|u;| 2 , where 
/i is the viscosity of the gel, and u) = — \e: {QQ T ) is the angular velocity of 
the fibrils with e being the permutation symbol. 

At steady temperature (an isothermal process, such as that in our experi- 
ments of Section 4.2) we have, from a Legendre transformation, e = ip + 9r), 
where rj is the entropy density of the system. In anticipation of the arguments 
to follow we now write the Hclmholtz free energy density as a sum of mechanical 
and chemical components, -0 = ip m + ip c where ip m = ip m (F, Q, X) is given by 
(21), which was written for a purely mechanical system. Therefore, the internal 
energy density rate satisfies 

e = V''m + i>c + Or). (27) 

The entropy density rate is governed by the Second Law of thermodynamics, 
which we write first as an equality: 

where the first term on the right hand-side is the entropy loss due to the heat 
sink, the second term is the entropy change due to the heat flux, and 7 is an 
entropy production term due to irreversible processes internal to the system. 
The statement of the Second Law as an inequality is 

7 > 0. (29) 

Multiplying (28) by 6, and combining it with (26) and (27), we have 
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Using ip m = ijj m (F, Q, X) as justified above, this can be expanded to 

The general hyperelastic constitutive law P = dtpm/dF reduces this inequality 
to 

-— : Q + Vc H ^ h 6*7 = 0. (30) 

Using (21) and (22) with ip(F,Q,X) replaced by $ m {F,Q,X) for the first 
term in (30), we have 



Nk B 8r 
AA 



Ar 

T 



{l-r/Lf 



4rg 
rL 



'0 



■F T F{Qe R ®Qe R ): 1 + -0, 



r(l-r /L) 2 



^0 

r 



+ 6» 7 = 0, (31) 



where M* = QMq and Mq = e R have been used. In a recent paper, Kuhl 
et al. (2005) have proposed the following first-order rate equation for the local 
evolution of fibril orientation: 



dM* 



1 



(M*-M max )M*-M max ] 



(32) 



dt t 

where r > is a relaxation time, and 7W max is the eigen vector corresponding 
to the maximum principal tensile stretch (see Remark 3). For the example in 
Section 4.3.1 we have JW max = e z . Equation (32) can be written in terms of Q 
by substituting M* = Qe R : 



Qe-R = -- [(e z -Qe R ) Qe R - e z ] . 



(33) 



Combining (33) and (31) we obtain 



Nk B 8r 
AA 



At 

T 



{l-r/Lf 



q-V9 



+ 6>7 = 0. 



F T F \ Qe R ® - [(e z -Qe R ) Qe R - e z ] 



ro 



rL r(l - r /LY 



ro 
r 



(34) 



Returning to the tensor product expansion for F, we use 1 = e R ®e R + e a ®e a + 
e z ®e z , and the fact that, for uniaxial tension, \ R = X a to write F T F = A^,l + 
(\ 2 Z — y? R )e Z <S> e z . Using this relation in (34) we get, after some manipulations, 



Nk B 8r 



AA 



At 
~L 



1 



(1-r/LY 



-1 - 



A 2 r 

■— (e Z -Qe R ) (e z -Qe R f-l 

T L 



4rg 
rL 



To 



r(l - r /L) 2 
q-V 



ro 
r 



6» 7 = 0. 



(35) 
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From (33) we have 



ez-Qen = — 

T 



(e z -Qe R ) 2 -l >0 (36) 



At t = the internal variable Q(t) = 1, leading to ez-Q(0)e R = . From this 
and (36) it follows that ez-Q{t)en > for all t > 0. Furthermore, for any locally 
remodelled orientation M* = Mz&z + MrBr, the unstretched and stretched 
lengths of the fibrils are in the ratio r^/r < 1 provided (Mz Xz) 2 + (MrA^) 2 > 1. 
Using these results the first term in (35) can be shown to be greater than zero 
provided the fibrils have rotated sufficiently for (MzAz) 2 + (Mr^r) 2 > 1 to 
hold, and r < L. 

This conclusion regarding the first term in (35) is not surprising. Its inter- 
pretation is that, for a given state of deformation F, the strain energy density 
at a point increases as Q locally rotates fibrils to align with the direction of 
maximum principal tensile stretch (see Arruda et al., 2005). In fact, this para- 
metric variation is a property of any fiber-reinforced material. The contributions 
of the remaining terms in (35) will now be studied towards understanding the 
dissipative character of biological remodelling. 

From the Second Law written as an inequality (29), it follows that the en- 
tropy production term is positive semi-definite. Assuming Fourier's Law of heat 
conduction, q = —K con V8, where K con is a positive semi-definite heat conduc- 
tion tensor, ensures that heat conduction results in a non-positive dissipation. 
Finally, since, as shown by our experiments (Section 4.2) re-orientation of fibrils 
happens by cell traction, the cells must consume their chemical free energy in 
this process: tp c < 0. Rewriting (35) with a reduced form for the first term on 
the left hand-side, 

-- 0. (37) 




Written in this form, the dissipation has several implications: 

1. Perhaps the most significant implication from the standpoint of develop- 
ment of models is that a purely mechanical theory is thcrmodynamically 
inadmissible for remodelling processes that stiffen the material. This is 
seen from (37) without the second and third terms, since the left hand- 
side is necessarily greater than zero if remodelling takes place. Either 
heat conduction or the chemical free energy changes must be considered 
to satisfy the dissipation relation. 

2. The chemical and mechanical action of cells needs to be considered in 
quantitative terms. While the chemical term is represented by ip c a con- 
stitutive model is still needed for it. Intra-cellular mechanical changes 
are also involved in the process (Section 4.2), and must be modelled by a 
separate free energy term. 
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3. The changes in chemical and mechanical free energy of the cells translates 
to entropy changes also. At this stage it is unclear whether the cells lower 
or increase their entropy by these processes. If the answer is a decrease, 
the overall positive entropy change must be sought further afield. 

4. While it is common to model biological tissues as isothermal, and even to 
ignore heat transport in them, (37) suggests that this is not necessarily 
allowable. Indeed it now appears that some heat flux must exist in the 
tissues for transport of the energetic and entropic byproducts of cellular 
activity. This term is indispensable if chcmcial free energy is ignored. 

6 Discussion and conclusion 

The hypothesis that biological systems attain stationary energy states with re- 
spect to changes in their microstructure is examined in this paper. A variational 
treatment can be applied to result in at least two, quite different, descriptions 
of remodelling: The first involving an evolution of the underlying material con- 
figuration, and the second described by internal variables. We note that the 
range of phenomena that can be described span from molecular to tissue scales. 

While the notion of stationary states of the free energy is examined, it is clear 
that dissipation must play a central role. This is highlighted by Section 5 and its 
result that stiffening of active biological materials requires both: consideration 
of chemical free energy and an entropic sink. This is a conclusion that merits 
deeper study. Attainment of true equilibrium states, in which all processes cease 
is dictated by dissipation. For this reason, and because biological systems have 
many mechanisms by which to dissipate energy, a very careful consideration of 
the Second Law and its consequences is essential to studies of remodelling. 

A Variational calculus for Case I remodelling 

Letting J* denote det(F*), ip denote the Helmholtz free energy density with re- 
spect to the reference configuration, and / denote the body force in the reference 
configuration, we rewrite the integrals in (6) over the reference configuration us- 
ing J*tp* = i> and J*f* = f. 
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Reverting to the remodelled configuration, 

4n[«*,«.] = /W^ + j^i 



de 



V de 
d« 



de 



- / / ■ -T^ dv * ~ t ■ — 7— dA 



de 



dfc E 



de 



on- 



Applying the chain rule and invoking stationarity with respect to variations in 

K, 



de 



TL[u*,k £ ] 



E=0 



- 1 f^.^ + j;^-^ 



\dK' de 



dF* ' de 



dV* 



E=0 



r -u: d £:^ + j: dr dK * 



dK' de 



dn 



f* ■ — -dV* + / t ■ —^-dA 



de 



6 dX* ■ de 



dn R 



dV* 



de 



c=0 



(38) 



e=0 



The derivatives with respect to e in (38) are obtained from (6): 



dK £ 
~de~ 



d5n 



E=0 



dX' V de 



dF! 



85k 



e=0 



<9X 



d/c, 



de 



5k. (39) 



r=0 



Substituting (39) in (38) gives 

d n[ u *,K £ ] 

e=0 



de 



, . r dip* _ 86k 
\dK ' 9F* ' 1 J 5JT 



5k dV* 



r -i (d£*_ dSn 

- 1 ( 7* dip* 85k dip 



dV* 



r dK : dX +J *dX* 



- I f* ■ SndV* - I t *■ SndA* 

n- an; 

Using 0J*/dK := d[det{K)]/dK = det(K)K- T , this reduces to 



de 



6=0 



K 



85k 

Ux 



dip* 



dF* 



(1-F*) 



8Sk 



dV* 



/ dip* dSn dip 



\dK' dX dX 



-5k dV* 



f* ■ SKdV* 



t* ■ SKdA* 



dn* t 
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Using standard manipulations of the scalar product of tensors, and the chain 
rule, 8(»)/8X* = {&{•)/ dX)K^ 1 , the preceding equation yields 



6=0 



85k f * t ^'* 85k . 

\8K : djr + dlc 



f* ■ 5ndV* - J t* ■ SkOA* 

n* an* 

Defining the configurational stress, £* := (dijj* /dK)K T , and introducing 
the first Piola-Kirchhoff stress, this can be written as 



de 



r=0 



9Sk , j, 88 n 85k . 

1> 1: T^ + ii-F ) P ■ ^ + — I dU* 



ax* 



d^* 

ax 5 



£«dU* -If- 8ndV* 



8X* 



t* ■ 8ndA* 



n* n* an* 

The Divergence Theorem allows this equation to be rewritten as 
d 



de 



e=0 



-Div* 1-F* P*+V*l + S 



8j'* 
8X~* 



SndV* 



(( 



1-F* 



+ S* ) • AT* - t 



8ndA* 



Enforcing equilibrium, ^nl/u*,/^] = 0, the localization principle and the 

£=0 

arbitrariness of Sk. give the following governing equations for configurational 
change: 



Div* P*+V*l + S* 



8JC 



- /* = in n* 
+ -0*1 + S*J AT* - r = on 90*. 
Using (3)i and (3)2 these equations are reducible to (7) and (8). 



(40) 
(41) 
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